function [primers]=primer_design
% defaults:
plength=21;
temp_range=[52 65];
gc_ratio=0.5;
overlap=300;
change_def=input('Change defaults?[0 1] ');
if change_def
    stop=0;
    while stop==0
        choice=input('Which parameter to change: length(0) temp_range(1) GCAT_ratio(2) overlap(3) ');
        switch choice
            case 0
                plength=input('Change length to: ');
            case 1
                temp_range=input('Change temp_range to: ');
            case 2
                gc_ratio=input('Change GCAT_ratio to: ');
            case 3
                overlap=input('Change overlap to: ');
        end
        stop=~input('Change more parameters?[0 1] ');
    end
end
min_read=500;
[file,pathid]=uigetfile('G:\Lab organization\MolBio\plasmids\maps\*.txt');

primers=struct;
primers.summary=struct;
primers.fwd=struct;
primers.rev=struct;

fid=fopen([pathid file]);
ini_seq=textscan(fid,'%s','delimiter','\n\r');
fclose(fid);
ini_seq=ini_seq{1};
seq={};
cnt=0;
for knd=2:size(ini_seq,1)
    cur_seq=ini_seq{knd};
    seq{1}(cnt+1:cnt+length(cur_seq))=cur_seq;
    cnt=cnt+length(cur_seq);
end
seq=seq{1};

seq_length=length(seq);
max_temp=max(temp_range);
min_temp=min(temp_range);
max_gc=gc_ratio+0.1;
min_gc=gc_ratio-0.1;

cnt=1;
fwd_primer_arr={};
for ind=1:length(seq)-plength
    wind=seq(ind:ind+(plength-1));
    g=findstr(wind,'G');
    c=findstr(wind,'C');
    a=findstr(wind,'A');
    t=findstr(wind,'T');
    obs_gc_ratio=(length(g)+length(c))/length(wind);
    last_nuc=wind(plength);
    if plength>14
        melt_temp=64.9+(41*((length(g)+length(c)-16.4)/plength));
    else
        melt_temp=(4*(length(g)+length(c)))+(2*(length(a)+length(t)));
    end
    run_g=findstr(wind,'GGG');
    run_c=findstr(wind,'CCC');
    run_a=findstr(wind,'AAA');
    run_t=findstr(wind,'TTT');
    run_g2=findstr(wind,'GGGG');
    run_c2=findstr(wind,'CCCC');
    run_a2=findstr(wind,'AAAA');
    run_t2=findstr(wind,'TTTT');
    if obs_gc_ratio>=min_gc & obs_gc_ratio<=max_gc
        if last_nuc=='G' | last_nuc=='C'
            if melt_temp>=min_temp & melt_temp<=max_temp
                if length(logical([run_g run_c run_a run_t run_g2 run_c2 run_a2 run_t2]))<2
                    fwd_primer_arr{cnt,1}=ind;
                    fwd_primer_arr{cnt,2}=wind;
                    fwd_primer_arr{cnt,3}=obs_gc_ratio;
                    fwd_primer_arr{cnt,4}=melt_temp;
                    cnt=cnt+1;
                end
            end
        end
    end
end

sequence_read=min_read+overlap;
min_n_prim=ceil(seq_length/min_read);

p_pos=[];
p_pos=cell2mat(fwd_primer_arr(:,1));
p_pos(end+1:end+length(p_pos))=p_pos+seq_length;
temp_arr={};
for ind=1:size(fwd_primer_arr,1)
    pos=p_pos(ind);
    count=0;
    temp_arr{ind,1}=pos;
    min_len=pos+min_read;
    max_len=pos+sequence_read;
    stop=0;
    while(stop==0)
        for gnd=1:(min_n_prim-1)
            if find(p_pos(p_pos>min_len & p_pos<max_len))
                tmp_pos=p_pos(p_pos>min_len & p_pos<max_len);
                pos_min=tmp_pos(1);
                pos_max=tmp_pos(end);
                min_len=pos_min+min_read;
                max_len=pos_max+sequence_read;
                for znd=1:length(tmp_pos)
                    if tmp_pos(znd)>seq_length
                        tmp_pos(znd)=tmp_pos(znd)-seq_length;
                    end
                end
                temp_arr{ind,gnd+1}=tmp_pos;
            end
            if gnd==(min_n_prim-1)
                stop=1;
            end
        end
    end
end


count2=0;
indi=[];
for ond=1:size(temp_arr,1)
    if ~isempty(temp_arr{ond,size(temp_arr,2)})
        count2=count2+1;
        indi(count2)=ond;
    end
end
fwd_choose_primer={};
fwd_choose_primer(1:length(indi),:)=temp_arr(indi,:);


for dnd=1:length(seq)
    if seq(length(seq)-(dnd-1))=='A'
        s2(1,dnd)='T';
    elseif seq(length(seq)-(dnd-1))=='T'
        s2(1,dnd)='A';
    elseif seq(length(seq)-(dnd-1))=='G'
        s2(1,dnd)='C';
    elseif seq(length(seq)-(dnd-1))=='C'
        s2(1,dnd)='G';
    end
end
cnt2=1;
rev_primer_arr={};
for jnd=1:length(s2)-plength
    wind=s2(jnd:jnd+(plength-1));
    g=findstr(wind,'G');
    c=findstr(wind,'C');
    a=findstr(wind,'A');
    t=findstr(wind,'T');
    obs_gc_ratio=(length(g)+length(c))/length(wind);
    last_nuc=wind(plength);
    if plength>14
        melt_temp=64.9+(41*((length(g)+length(c)-16.4)/plength));
    else
        melt_temp=(4*(length(g)+length(c)))+(2*(length(a)+length(t)));
    end
    run_g=findstr(wind,'GGG');
    run_c=findstr(wind,'CCC');
    run_a=findstr(wind,'AAA');
    run_t=findstr(wind,'TTT');
    run_g2=findstr(wind,'GGGG');
    run_c2=findstr(wind,'CCCC');
    run_a2=findstr(wind,'AAAA');
    run_t2=findstr(wind,'TTTT');
    if obs_gc_ratio>=0.9 & obs_gc_ratio<=1.1
        if last_nuc=='G' | last_nuc=='C'
            if melt_temp>=55 & melt_temp<=70
                if length(logical([run_g run_c run_a run_t run_g2 run_c2 run_a2 run_t2]))<2
                    rev_primer_arr{cnt2,1}=jnd;
                    rev_primer_arr{cnt2,2}=wind;
                    rev_primer_arr{cnt2,3}=obs_gc_ratio;
                    rev_primer_arr{cnt2,4}=melt_temp;
                    cnt2=cnt2+1;
                end
            end
        end
    end
end
p_pos=[];
p_pos=cell2mat(rev_primer_arr(:,1));
p_pos(end+1:end+length(p_pos))=p_pos+seq_length;
temp_arr={};
for ind=1:size(rev_primer_arr,1)
    pos=p_pos(ind);
    count=0;
    temp_arr{ind,1}=pos;
    min_len=pos+min_read;
    max_len=pos+sequence_read;
    stop=0;
    while(stop==0)
        for gnd=1:(min_n_prim-1)
            if find(p_pos(p_pos>min_len & p_pos<max_len))
                tmp_pos=p_pos(p_pos>min_len & p_pos<max_len);
                pos_min=tmp_pos(1);
                pos_max=tmp_pos(end);
                min_len=pos_min+min_read;
                max_len=pos_max+sequence_read;
                for znd=1:length(tmp_pos)
                    if tmp_pos(znd)>seq_length
                        tmp_pos(znd)=tmp_pos(znd)-seq_length;
                    end
                end
                temp_arr{ind,gnd+1}=tmp_pos;
            end
            if gnd==(min_n_prim-1)
                stop=1;
            end
        end
    end
end

count2=0;
indi=[];
for ond=1:size(temp_arr,1)
    if ~isempty(temp_arr{ond,size(temp_arr,2)})
        count2=count2+1;
        indi(count2)=ond;
    end
end
rev_choose_primer={};
rev_choose_primer(1:length(indi),:)=temp_arr(indi,:);

primers.summary.seq_length=seq_length;
primers.summary.max_temp=max_temp;
primers.summary.min_temp=min_temp;
primers.summary.max_GC=max_gc;
primers.summary.min_GC=min_gc;
primers.summary.overlap=overlap;
primers.summary.min_n_prim=min_n_prim;
primers.fwd.full_primer_arr=fwd_primer_arr;
primers.fwd.overlapping_primers=fwd_choose_primer;
primers.rev.full_primer_arr=rev_primer_arr;
primers.rev.overlapping_primers=rev_choose_primer;

for fnd=1:size(fwd_choose_primer,1)
    for tnd=1:size(fwd_choose_primer,2)
        best_sel_fwd(fnd,tnd)=fwd_choose_primer{fnd,tnd}(end);
    end
end
for fnd=1:size(rev_choose_primer,1)
    for tnd=1:size(rev_choose_primer,2)
        best_sel_rev(fnd,tnd)=rev_choose_primer{fnd,tnd}(end);
    end
end

primers.fwd.best_sel=best_sel_fwd;
primers.rev.best_sel=best_sel_rev;

if size(fwd_choose_primer,2)>size(rev_choose_primer,2)
    best_set_len_fwd=sum(abs(diff(primers.fwd.best_sel')));
    best_set_ind=find(max(best_set_len_fwd));
    primers.best_set_fwd.loc =primers.fwd.best_sel(best_set_ind,:);
    list_loc=primers.fwd.full_primer_arr(:,1);
    seq=char();
    for mnd=1:length(primers.best_set_fwd.loc)
        curr_prim_loc=primers.best_set_fwd.loc(mnd);
        for nnd=1:length(list_loc)
            curr_loc=list_loc{nnd};
            if curr_loc==curr_prim_loc
                seq(mnd,:)=char(primers.fwd.full_primer_arr(nnd,2));
            end
        end
    end
    primers.best_set_fwd.seq =seq;
elseif size(fwd_choose_primer,2)<size(rev_choose_primer,2)
    best_set_len_rev=sum(abs(diff(primers.rev.best_sel')));
    best_set_ind=find(max(best_set_len_rev));
    primers.best_set_rev.loc =primers.rev.best_sel(best_set_ind,:);
    list_loc=primers.rev.full_primer_arr(:,1);
    seq=char();
    for mnd=1:length(primers.best_set_rev.loc)
        curr_prim_loc=primers.best_set_rev.loc(mnd);
        for nnd=1:length(list_loc)
            curr_loc=list_loc{nnd};
            if curr_loc==curr_prim_loc
                seq(mnd,:)=char(primers.rev.full_primer_arr(nnd,2));
            end
        end
    end
    primers.best_set_rev.seq =seq;
elseif size(fwd_choose_primer,2)==size(rev_choose_primer,2)
    best_set_len_fwd=sum(abs(diff(primers.fwd.best_sel')));
    best_set_ind_fwd=find(max(best_set_len_fwd));
    primers.best_set_fwd.loc =primers.fwd.best_sel(best_set_ind_fwd,:);
    
    best_set_len_rev=sum(abs(diff(primers.rev.best_sel')));
    best_set_ind_rev=find(max(best_set_len_rev));
    primers.best_set_rev.loc =primers.rev.best_sel(best_set_ind_rev,:);
    
    list_loc=primers.fwd.full_primer_arr(:,1);
    seq=char();
    for mnd=1:length(primers.best_set_fwd.loc)
        curr_prim_loc=primers.best_set_fwd.loc(mnd);
        for nnd=1:length(list_loc)
            curr_loc=list_loc{nnd};
            if curr_loc==curr_prim_loc
                seq(mnd,:)=char(primers.fwd.full_primer_arr(nnd,2));
            end
        end
    end
    primers.best_set_fwd.seq =seq;
    list_loc=primers.rev.full_primer_arr(:,1);
    seq=char();
    for mnd=1:length(primers.best_set_rev.loc)
        curr_prim_loc=primers.best_set_rev.loc(mnd);
        for nnd=1:length(list_loc)
            curr_loc=list_loc{nnd};
            if curr_loc==curr_prim_loc
                seq(mnd,:)=char(primers.rev.full_primer_arr(nnd,2));
            end
        end
    end
    primers.best_set_rev.seq =seq;
end
